home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / b1nuI.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-11-24  |  7.8 KB  |  316 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: b1nuI.c,v 1.4 85/08/22 16:51:13 timo Exp $
  5. */
  6.  
  7. /* Multi-precision integer arithmetic */
  8.  
  9. #include "b.h"
  10. #include "b1obj.h"
  11. #include "b1num.h"
  12. #include "b0con.h"
  13. #include "b3err.h"
  14.  
  15. /*
  16.  * Number representation:
  17.  * ======================
  18.  *
  19.  * (Think of BASE = 10 for ordinary decimal notation.)
  20.  * A number is a sequence of N "digits" b1, b2, ..., bN
  21.  * where each bi is in {0..BASE-1}, except for negative numbers,
  22.  * where bN = -1.
  23.  * The number represented by b1, ..., bN is
  24.  *      b1*BASE**(N-1) + b2*BASE**(N-2) + ... + bN .
  25.  * The base BASE is chosen so that multiplication of two positive
  26.  * integers up to BASE-1 can be multiplied exactly using double
  27.  * precision floating point arithmetic.
  28.  * Also it must be possible to add two long integers between
  29.  * -BASE and +BASE (exclusive), giving a result between -2BASE and
  30.  * +2BASE.
  31.  * BASE must be even (so we can easily decide whether the whole
  32.  * number is even), and positive (to avoid all kinds of other trouble).
  33.  * Presently, it is restricted to a power of 10 by the I/O-conversion
  34.  * routines (file "b1nuC.c").
  35.  *
  36.  * Canonical representation:
  37.  * bN is never zero (for the number zero itself, N is zero).
  38.  * If bN is -1, b[N-1] is never BASE-1 .
  39.  * All operands are assumed te be in canonical representation.
  40.  * Routine "int_canon" brings a number in canonical representation.
  41.  *
  42.  * Mapping to C objects:
  43.  * A "digit" is an integer of type "digit", probably an "int".
  44.  * A number is represented as a "B-integer", i.e. something
  45.  * of type "integer" (which is actually a pointer to some struct).
  46.  * The number of digits N is extracted through the macro Length(v).
  47.  * The i-th digit is extracted through the macro Digit(v,N-i).
  48.  * (So in C, we count in a backwards direction from 0 ... n-1 !)
  49.  * A number is created through a call to grab_num(N), which sets
  50.  * N zero digits (thus not in canonical form!).
  51.  */
  52.  
  53.  
  54. /*
  55.  * Bring an integer into canonical form.
  56.  * Make a SmallInt if at all possible.
  57.  * NB: Work done by int_canon is duplicated by mk_integer for optimization;
  58.  *     if the strategy here changes, look at mk_integer, too!
  59.  */
  60.  
  61. Visible integer int_canon(v) integer v; {
  62.     register int i;
  63.  
  64.     if (IsSmallInt(v)) return v;
  65.  
  66.     for (i = Length(v) - 1; i >= 0 && Digit(v,i) == 0; --i)
  67.         ;
  68.  
  69.     if (i < 0) {
  70.         release((value) v);
  71.         return int_0;
  72.     }
  73.  
  74.     if (i == 0) {
  75.         digit dig = Digit(v,0);
  76.         release((value) v);
  77.         return (integer) MkSmallInt(dig);
  78.     }
  79.  
  80.     if (i > 0 && Digit(v,i) == -1) {
  81.         while (i > 0 && Digit(v, i-1) == BASE-1) --i;
  82.         if (i == 0) {
  83.             release((value) v);
  84.             return (integer) MkSmallInt(-1);
  85.         }
  86.         if (i == 1) {
  87.             digit dig = Digit(v,0) - BASE;
  88.             release((value) v);
  89.             return (integer) MkSmallInt(dig);
  90.         }
  91.         Digit(v,i) = -1;
  92.     }
  93.  
  94.     if (i+1 < Length(v)) return (integer) regrab_num((value) v, i+1);
  95.  
  96.     return v;
  97. }
  98.  
  99.  
  100. /* General add/subtract subroutine */
  101.  
  102. typedef double twodigit; /* Might be long on 16 bit machines */
  103.     /* Should be in b0con.h */
  104.  
  105. Hidden twodigit fmodulo(x, y) twodigit x, y; {
  106.     return x - y * (twodigit) floor((double)x / (double)y);
  107. }
  108.  
  109. Visible Procedure dig_gadd(to, nto, from, nfrom, ffactor)
  110.     digit *to, *from; intlet nto, nfrom; digit ffactor; {
  111.     twodigit carry= 0;
  112.     twodigit factor= ffactor;
  113.     digit save;
  114.  
  115.     nto -= nfrom;
  116.     if (nto < 0)
  117.         syserr(MESS(1000, "dig_gadd: nto < nfrom"));
  118.     for (; nfrom > 0; ++to, ++from, --nfrom) {
  119.         carry += *to + *from * factor;
  120.         *to= save= fmodulo(carry, (twodigit)BASE);
  121.         carry= (carry-save) / BASE;
  122.     }
  123.     for (; nto > 0; ++to, --nto) {
  124.         if (carry == 0)
  125.             return;
  126.         carry += *to;
  127.         *to= save= fmodulo(carry, (twodigit)BASE);
  128.         carry= (carry-save) / BASE;
  129.     }
  130.     if (carry != 0)
  131.         to[-1] += carry*BASE; /* Assume it's -1 */
  132. }
  133.  
  134.  
  135. /* Sum or difference of two integers */
  136. /* Should have its own version of dig-gadd without double precision */
  137.  
  138. Visible integer int_gadd(v, w, factor) integer v, w; intlet factor; {
  139.     struct integer vv, ww;
  140.     integer s;
  141.     int len, lenv, i;
  142.  
  143.     FreezeSmallInt(v, vv);
  144.     FreezeSmallInt(w, ww);
  145.     lenv= len= Length(v);
  146.     if (Length(w) > len)
  147.         len= Length(w);
  148.     ++len;
  149.     s= (integer) grab_num(len);
  150.     for (i= 0; i < lenv; ++i)
  151.         Digit(s, i)= Digit(v, i);
  152.     for (; i < len; ++i)
  153.         Digit(s, i)= 0;
  154.     dig_gadd(&Digit(s, 0), len, &Digit(w, 0), Length(w), (digit)factor);
  155.     return int_canon(s);
  156. }
  157.  
  158.  
  159. /* Product of two integers */
  160.  
  161. Visible integer int_prod(v, w) integer v, w; {
  162.     int i;
  163.     integer a;
  164.     struct integer vv, ww;
  165.  
  166.     if (v == int_0 || w == int_0) return int_0;
  167.     if (v == int_1) return (integer) Copy(w);
  168.     if (w == int_1) return (integer) Copy(v);
  169.  
  170.     FreezeSmallInt(v, vv);
  171.     FreezeSmallInt(w, ww);
  172.  
  173.     a = (integer) grab_num(Length(v) + Length(w));
  174.  
  175.     for (i= Length(a)-1; i >= 0; --i)
  176.         Digit(a, i)= 0;
  177.     for (i = 0; i < Length(v) && !interrupted; ++i)
  178.         dig_gadd(&Digit(a, i), Length(w)+1, &Digit(w, 0), Length(w), 
  179.             Digit(v, i));
  180.  
  181.     return int_canon(a);
  182. }
  183.  
  184.  
  185. /* Compare two integers */
  186.  
  187. Visible relation int_comp(v, w) integer v, w; {
  188.     int sv, sw;
  189.     register int i;
  190.     struct integer vv, ww;
  191.  
  192.     /* 1. Compare pointers and equal SmallInts */
  193.     if (v == w) return 0;
  194.  
  195.     /* 1a. Handle SmallInts */
  196.     if (IsSmallInt(v) && IsSmallInt(w))
  197.         return SmallIntVal(v) - SmallIntVal(w);
  198.     FreezeSmallInt(v, vv);
  199.     FreezeSmallInt(w, ww);
  200.  
  201.     /* 2. Extract signs */
  202.     sv = Length(v)==0 ? 0 : Digit(v,Length(v)-1)<0 ? -1 : 1;
  203.     sw = Length(w)==0 ? 0 : Digit(w,Length(w)-1)<0 ? -1 : 1;
  204.  
  205.     /* 3. Compare signs */
  206.     if (sv != sw) return (sv>sw) - (sv<sw);
  207.  
  208.     /* 4. Compare sizes */
  209.     if (Length(v) != Length(w))
  210.         return sv * ( (Length(v)>Length(w)) - (Length(v)<Length(w)) );
  211.  
  212.     /* 5. Compare individual digits */
  213.     for (i = Length(v)-1; i >= 0 && Digit(v,i) == Digit(w,i); --i)
  214.         ;
  215.  
  216.     /* 6. All digits equal? */
  217.     if (i < 0) return 0;  /* Yes */
  218.  
  219.     /* 7. Compare leftmost different digits */
  220.     if (Digit(v,i) < Digit(w,i)) return -1;
  221.  
  222.     return 1;
  223. }
  224.  
  225.  
  226. /* Construct an integer out of a floating point number */
  227.  
  228. #define GRAN 8    /* Granularity used when requesting more storage */
  229.         /* MOVE TO MEM! */
  230. Visible integer mk_int(x) double x; {
  231.     register integer a;
  232.     integer b;
  233.     register int i, j;
  234.     int negate;
  235.  
  236.     if (MinSmallInt <= x && x <= MaxSmallInt)
  237.         return (integer) MkSmallInt((int)x);
  238.  
  239.     a = (integer) grab_num(1);
  240.     negate = x < 0 ? 1 : 0;
  241.     if (negate) x = -x;
  242.  
  243.     for (i = 0; x != 0; ++i) {
  244.         double z = floor(x/BASE);
  245.         digit save = Modulo((digit)(x-z*BASE), BASE);
  246.         if (i >= Length(a)) {
  247.             a = (integer) regrab_num((value) a, Length(a)+GRAN);
  248.             for (j = Length(a)-1; j > i; --j)
  249.                 Digit(a,j) = 0;    /* clear higher digits */
  250.         }
  251.         Digit(a,i) = save;
  252.         x = floor((x-save)/BASE);
  253.     }
  254.  
  255.     if (negate) {
  256.         b = int_neg(a);
  257.         release((value) a);
  258.         return b;
  259.     }
  260.  
  261.     return int_canon(a);
  262. }
  263.  
  264. /* Construct an integer out of a C int.  Like mk_int, but optimized. */
  265.  
  266. Visible value mk_integer(x) int x; {
  267.     if (MinSmallInt <= x && x <= MaxSmallInt) return MkSmallInt(x);
  268.     return (value) mk_int((double)x);
  269. }
  270.  
  271.  
  272. /* Efficiently compute 10**n as a B integer, where n is a C int >= 0 */
  273.  
  274. Visible integer int_tento(n) int n; {
  275.     integer i;
  276.     digit msd = 1;
  277.     if (n < 0) syserr(MESS(1001, "int_tento(-n)"));
  278.     if (n < tenlogBASE) {
  279.         while (n != 0) msd *= 10, --n;
  280.         return (integer) MkSmallInt(msd);
  281.     }
  282.     i = (integer) grab_num(1 + (int)(n/tenlogBASE));
  283.     n %= tenlogBASE;
  284.     while (n != 0) msd *= 10, --n;
  285.     Digit(i, Length(i)-1) = msd;
  286.     return i;
  287. }
  288.  
  289. #ifdef NOT_USED
  290. /* Approximate ceiling(10 log abs(u/v)), as C int.
  291.    It only works for v > 0, u, v both integers.
  292.    The result may be one too large or too small */
  293.  
  294. Visible int scale(u, v) integer u, v; {
  295.     int s;
  296.     double z;
  297.     struct integer uu, vv;
  298.  
  299.     if (Msd(v) <= 0) syserr(MESS(1002, "scale(u,v<=0)"));
  300.     if (u == int_0) return 0; /* `Don't care' case */
  301.     FreezeSmallInt(u, uu);
  302.     FreezeSmallInt(v, vv);
  303.     s = (Length(u) - Length(v)) * tenlogBASE;
  304.     if (Digit(u, Length(u)-1) >= 0) z = Digit(u, Length(u)-1);
  305.     else {
  306.         s -= tenlogBASE;
  307.         if (Length(u) == 1) z = 1;
  308.         else z = BASE - Digit(u, Length(u)-2);
  309.     }
  310.     z /= Digit(v, Length(v)-1);
  311.     while (z >= 10) z /= 10, ++s;
  312.     while (z < 1) z *= 10, --s;
  313.     return s;
  314. }
  315. #endif NOT_USED
  316.